
delta_1 = 1;
ind_1   = find(delta_grid == delta_1);

s_hat  = fn_s_hat(R1,lambda,s_min,s_max,phi); % Computes s_hat, which is indepdent of delta
dq_Fdd = fn_dq(delta_grid,R1,D0,lambda,prob,dist_s,s_min,s_max,dist_d,d_min,d_max,phi,delta_grid_stp); % Computes partial of q wrt delta

s_star_1       = fn_s_star(delta_1,R1,D0,lambda,s_min,s_max,dist_d,d_min,d_max,phi);
s_star_0       = fn_s_star(0,R1,D0,lambda,s_min,s_max,dist_d,d_min,d_max,phi);
q_F_0          = dist_s.cdf(s_hat) + prob*(dist_s.cdf(s_star_0) - dist_s.cdf(s_hat));
q_F_1          = dist_s.cdf(s_hat) + prob*(dist_s.cdf(s_star_1) - dist_s.cdf(s_hat));
q_FF_1         = dist_s.cdf(s_hat);

dq_Fdd_1                  = dq_Fdd(ind_1);
frac_dq_Fdd_1             = dq_Fdd_1/q_F_1;
[q_T_1, q_T_F_1, E_tax_F_1, E_tax_T_F_1, E_MC_tax_F_1, E_MC_tax_T_F_1, max_MC_tax_1] = ...
    fn_E_tax(R1,delta_1,D0,chi,phi,s_min,s_max,dist_s,dist_d,d_min,d_max,prob,s_hat,s_star_1,kappa);

[ins_dep_1, unins_dep_1] = fn_insured_deposits(delta_1,R1,dist_d,d_min,d_max);
[ins_acc_1, unins_acc_1] = fn_insured_accounts(delta_1,R1,dist_d,d_min,d_max);

T_star_1    = fn_tax(delta_1,R1,D0,chi,s_star_1,phi,dist_d,d_min,d_max); % T(s^*)

f_Cdif_i    = @(D0i) fn_Cdif_i(D0i,R1,delta_1,D0,chi,lambda,y1,y2,phi,dist_d,d_min,d_max,s_star_1).*dist_d.pdf(D0i);
Cdif_dep_1  = integral(f_Cdif_i,d_min,d_max);
Cdif_tax_1  = fn_kappa_gross(T_star_1,kappa);
Cdif_1      = Cdif_dep_1 + Cdif_tax_1;
